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Abstract 

In this paper, we propose an efficient and spectrally accurate numerical 
method for computing the dynamics of rotating Bose-Einstein condensates 
(BEC) in two dimensions (2D) and 3D based on the Gross-Pitaevskii equa- 
tion (GPE) with an angular momentum rotation term. By applying a time- 
splitting technique for decoupling the nonlinearity and properly using the 
alternating direction implicit (ADI) technique for the coupling in the angu- 
lar momentum rotation term in the GPE, at every time step, the GPE in 
rotational frame is decoupled into a nonlinear ordinary differential equation 
(ODE) and two partial differential equations with constant coefficients. This 
allows us to develop new time-splitting spectral (TSSP) methods for comput- 
ing the dynamics of BEC in a rotational frame. The new numerical method is 
explicit, unconditionally stable, and of spectral accuracy in space and second 
order accuracy in time. Moreover, it is time reversible and time transverse 
invariant, and conserves the position density in the discretized level if the 
GPE does. Extensive numerical results are presented to confirm the above 
properties of the new numerical method for rotating BEC in 2D & 3D. 

Key Words: Rotating Bose-Einstein condensates, Gross-Pitaevskii equation, an- 
gular momentum rotation, time spitting. 



1 Introduction 

Since the first experimental creation of a quantized vortex in a gaseous Bose-Einstein 
condensate (BEC) [29, 1, 30, 35], there has been significantly experimental and 
theoretical advances in the field of research [3, 34, 2, 11, 5, 13, 23, 14, 17, 19, 22]. 
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Several experimental methods of vortex creation are currently in use, including 
phase imprinting [30, 40], cooling of a rotating normal gas [24], and conversion of 
spin angular momentum into orbital angular momentum by reversal of the magnetic 
bias field in an loffc-Pritchard trap [26, 27, 33]. The topic of this paper is to propose 
an efficient and spectrally accurate numerical method for studying quantized vortex 
dynamics in a BEC by imposing a laser beam rotating with an angular velocity on 
the magnetic trap holding the atoms to create a harmonic anisotropic potential. 

At temperatures T much smaller than the critical condensation temperature Tc, 
under mean field theory, the properties of a BEC in a rotational frame are modelled 
by the well-known time-dependent Gross-Pitaevskii equation (GPE) with an angular 
momentum rotation term [3, 14, 22]: 

av^^ ^ (-^^' + + ^^ol^r - ^L^j ^(x, t), X e M^ t > 0, 

(1.1) 

where x = (x, zY is the Cartesian coordinate vector, ■?/'(x, t) is the complex-valued 
macroscopic wave function, m is the atomic mass, h is the Planck constant, N is the 
number of atoms in the condensate, Vt is the angular velocity of rotating laser beam, 
\^(x) = y i^x^"^ + '^yV^ + '^1^'^^ with ojx, u!y and cUz being the trap frequencies in x-, 
y- and ^-direction respectively. Uq = describes the interaction between atoms 

in the condensate with the s-wave scattering length, and Lz is the z-component of 
the angular momentum. It is convenient to normalize the wave function by requiring 

m;t)r--= [m^,t)\'d^^i. (1.2) 

After proper nondimensionalization and dimension reduction, we can obtain the 
following dimensionless GPE with an angular momentum rotation term in the d- 
dimensions {d — 2,3) [13, 5, 41]: 

i^^^^^^-lv'^l^ + Va{^)ij + PM'i^-^L,^: xeR^ t>0, (1.3) 
^/'(x,0) = Vo(x), xeM^ with Uof-^f |^o(x)|'(ix= 1, (1.4) 

where — —i{xdy — ydx) and 

, / + /2, d = 2, 

^d(x) = M 2 2 ^ 2 2^ 2 2\ /o , ^ (1-5) 

with 7^., and 7^ being constants. 

In order to study effectively the dynamics of BEC, especially in the strong re- 
pulsive interaction regime, i.e. /5d ^ 1 in (1.3), an efficient and accurate numerical 
method is one of the key issues. For non-rotating BEC, i.e. Q = in (1.3), many 
efficient and spectrally accurate numerical methods were proposed in the literatures 
[7, 6, 11, 4, 12], and they were demonstrated that they are much better than the 
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low-order finite difference metliods [15, 32, 31, 16]. Tlius tliey were applied to study 
collapse and explosion of BEC in 3D [8] and multi-component BEC [4] which are the 
very challenging problems in numerical simulation of BEC. Due to the appearance 
of the angular momentum rotation term in the GPE (1.3), new numerical difficulties 
are introduced. Currently, the numerical methods used in the physics literature for 
studying dynamics of rotating BEC remain limited [25, 39], and they usually are 
low-order finite difference methods. Recently, some efficient and accurate numerical 
methods were designed for computing dynamics of rotating BEC. For example, Bao, 
Du and Zhang [5] proposed a numerical method for computing dynamics of rotating 
BEC by applying a time-splitting technique for decoupling the nonlinearity in the 
GPE and adopting the polar coordinates or cylindrical coordinates so as to make 
the coefficient of the angular momentum rotation term constant. The method is 
time reversible, unconditionally stable, implicit but can be solved very efficiently, 
and conserves the total density. It is of spectral accuracy in transverse direction, but 
usually of second or fourth-order accuracy in radial direction. Another numerical 
method is the leap-frog spectral method used for studying vortex lattice dynamics 
in rotating BEC [41]. This method is exphcit, time reversible, of spectral accuracy 
in space and second order accuracy in time. But it has a stability constraint for 
time step [41]. The aim of this paper is to develop a numerical method which enjoys 
advantages of both the above two numerical methods. That is to say, the method 
is explicit, unconditionally stable, time reversible, time transverse invariant, and of 
spectral accuracy in space. We shall present such an efficient, unconditionally stable 
and accurate numerical method for discretizing the GPE in a rotational frame by 
applying a time-splitting technique and an ADI technique, and constructing appro- 
priately spectral basis functions. 

The paper is organized as follows. In section 2, wc review some properties of 
GPE in a rotational frame (1.3) including conservation laws and analytical solutions 
of condensate widths. In section 3, we present a new time-splitting Fourier pseu- 
dospectral method for efficient and accurate simulation of GPE (1.3) in 2D & 3D. In 
section 4, extensive numerical results are reported to demonstrate the efficiency and 
spectral resolution in space of our new numerical method. Finally some conclusions 
are drawn in section 5. 

2 Some properties of the GPE 

For the convenience of the reader, in this section, we will review some properties 
of the GPE with an angular momentum rotation term (1.3), which will be used to 
test our new numerical method proposed in the next section. First of all, the GPE 
(1.3) is time reversible and time transverse invariant. Second, it has at least two 
important invariants which are the normalization of the wave function 

Ni^l^)^ [ \^{^,t)\'d^= [ |V'(x,0)|2dx = Ar(V'o) = l, t>0, (2.1) 
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and the energy 



= S/3„f^(^o), t>0, (2.2) 

where /* and Ini(/) denote the conjugate and the imaginary part of the function 
/ respectively. Third, it was proven that at least for radial symmetric trap in 2D 
or cylindrical symmetric trap in 3D, i.e., — ly in (1-3), the angular momentum 
expectation and energy for non- rotating part are conserved [5], that is, for any given 
initial data 'V'o(x) in (1.4), 

(L,)(i) = (L,)(0), £;^,o(V) = £^/3,o(V'o), i>0, (2.3) 

where the angular momentum expectation which is a measure of the vortex flux is 

defined as 

{L^){t):^[ i/j*{yi,t)L^ijj{x,t) dx^i [ 'ijj*U,t)(ydx - xdyU(x,t)d^, t>0. 

(2.4) 

Other very useful quantities characterizing the dynamics of rotating BEG in 2D are 

the condensate width defined as 

5x{t)=[ a:2|V^(x,t)|2dx, 5y{t)=[ y2|^(x,t)|2(ix, 5,(t) = 5,(t) + 5,(t). (2.5) 

As proven in [5], in 2D with a radial symmetric trap, i.e., d = 2 and — ly '■— Ir 
in (1.3), for any initial data ipQ^x^y) in (1.4), we have for any t > 

= EpM + ^{Lm _ ^^^(27.t)] + 5f cos(27.t) + sin(27.t), (2.6) 

7r 27. 

where df^ := 5.(0) = 5^(0) + 5y{Q) and 5^^^ := (^.(0) = (5^(0) + 4(0) with for a = x 
or y 

'^a(0) = <5f = / a^|V'o(x)|'dx, 

^a(O) = 6^^^ = 2 [ a [-filV^oP {xdy - yd.,) a + lm (rodai^o)] d^. 

Furthermore, when the initial data ilJo{x,y) in (1.4) satisfies 

^jjo{x,y) ^ f{r)e''^^ with meZ and /(O) = when ruy^O, (2.7) 
with (r, 9) being the polar coordinates in 2D, we have, for any t > 0, 

= ^A^(^o) + [1 - eos(27.t)] + cos(27.t) + ^ sin(27.t). (2.8) 

These immediately imply that 5.(t) is a periodic function with angular frequency 
doubling the trapping frequency in 2D with a radial symmetric trap, and also Sx{t) 
and Sy{t) are periodic functions with frequency doubling the trapping frequency 
provided that the initial data satisfies (2.7). 
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3 A time-splitting pseudospectral method for ro- 
tating BEC 

In this section, we will present an explicit, unconditionally stable and spectrally 
accurate numerical method to solve the GPE (1.3) for dynamics of rotating BEC. 

Due to the external trapping potential Ki(x), the solution ■(/'(x, t) of (1.3)-(1.4) 
decays to zero exponentially fast when |x| oo. Thus in practical computation, we 
always truncate the problem (1.3)-(1.4) into a bounded computational domain with 
homogeneous Dirichlet boundary condition: 

i 9tV'(x, t) = V'V + [Vi^) - iW(x.)] ip + PaW'jjj - QL,V, x e fi,, t > 0, (3.1) 
V'(x, t)^0, X e r = dn^, t > 0, (3.2) 
V'(x,0) =V'o(x), xeQ,; (3.3) 

where W^(x) > corresponds to a localized loss term [36] and V{yi) = V^(x) + Vp{x) 
with Vp{'x) > a conservative repulsive pinning potential [36]. Here we choose 
= [a, &] X [c,d] in 2D, and resp., Qx = [o, &] x [c,d\ x [e, /] in 3D, with |a|, b, 
|c|, d, \e\ and / sufficiently large. The use of more sophisticated radiation boundary 
conditions is an interesting topic that remains to be examined in the future. 

3.1 Time-splitting 

Wc choose a time step size At > 0. For n = 0, 1, 2, ■ ■ -, from time t = tn = nAt to 

t = tn+i = tn + At, the GPE (3.1) is first solved in two splitting steps. One solves 
first 

t dM^, t) = V2^(x, t) - l^L,V^(x, t) (3.4) 
for the time step of length At, followed by solving 

i dM^,t) = [V^(x) -^iy(x)]^(x,t) +/5,|^(x,t)|2^(x,t), (3.5) 

for the same time step. For t e [tn,tn+i], multiplying (3.5) by the conjugate of 

ip, we get 

i ri^,t)dM^,t) = [y(x) - iW{^)] |V'(x,t)|^ + /3<i|V'(x,t)|^ (3.6) 
Subtracting the conjugate of (3.6) from (3.6) and multiplying by —i, one obtains 

^|^(x,t)|2 = rdti^ + i^dtr = -2iy(x)|^(x,t)|2. (3.7) 

Solving (3.7), we get 

|^(x,t)|2 = e-^^W(*-*")|^(x,y|2, < t < t„+i. (3.8) 
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Substituting (3.8) into (3.5), we obtain 

i dM^,t) = [l-(x) -zVF(x)]^(x,t) +/3,e-2^W(*-*")|V^(x,y|2V;(x,t). (3.9) 
Integrate (3.9) from tn to t, we find for x e and tn <t < tn+i- 

' g-i[y(x)+/3d|v-(x,t„)|^]{t-tn) ^(x,i„), W{^) = 0, 

V'(x,i) = < 

(3.10) 

To discretize (3.4) in space, compared with non-rotating BEC [7, 11, 6], i.e. 1] = 
in (1.3), the main difficulty is that the coefficients in are not constants which 
cause big trouble in applying sine or Fourier pseudospectral discretization. Due 
to the special structure in the angular momentum rotation term (1.5), wc will 
apply the alternating direction implicit (ADl) technique and decouple the operator 
— — Q.Lz into two one dimensional operators such that each operator becomes 
a summation of terms with constant coefficients in that dimension. Therefore, they 
can be discrctizcd in space by Fourier pseudospectral method and the ODEs in phase 
space can be integrated analytically. The details for discretizing (3.4) in 2D & 3D 
will be presented in the next two subsections respectively. 

3.2 Discretization in 2D 

When c? = 2 in (3.4), we choose mesh sizes A,x > and Ay > with Ax = {h — a)/M 
and Ay = [d — c)/N for M and N even positive integers, and let the grid points be 

Xj^a + jAx, J =0,1,2,---,M; yk^c + kAy, k ^ 0,1,2, N. 

Let ip'jj^ be the approximation of ip{xj,yk,tn) and ■0" be the solution vector with 
component tpji^. 

From time t = tn to t = tn+i, we solve (3.4) first 

i dtip{x, t) = -^9^^V(x, t) - iQyd:,ijj{yi, t), (3.11) 

for the time step of length At, followed by solving 

i dtip{-x, t) — —^dyyip{-x., t) + iflxdyip{-K, t), (3-12) 

for the same time step. The detailed discretizations of (3.11) and (3.12) are shown 
in Appendix A. 



3.3 Discretization in 3D 

When d = 3 in (3.4), we choose mesh sizes Ax > 0, Ay > and Az > with 
Ax = {b-a)/M, Ay = {d-c)/N and Az = {f-e)/L for M, N and L even positive 
integers, and let the grid points be 

xj jAx, <j < M; yk^c + kAy, < k < N; zi^e + lAz, < I < L. 
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Let ipj'f.i be the approximation of il'{xj, yk, zi, tn) and ip"' be the solution vector with 
component ipji^i- 

Similar as those for 2D case, from time t — tn to t — tn+i, we solve (3.4) first 

i 9tV'(x, t) = (-^9^x - ^d^z - i^ydx^ V'(x, t), (3.13) 



for the time step of length At, followed by solving 

1 . 1 



i dtipi^, t) = (^-^dyy - -dzz + iQxdy^ ^(x, t), 



(3.14) 



for the same time step. The detailed discretizations of (3.13) and (3.14) are shown 
in Appendix B. 



3.4 Stability 

We define the usual discrete Z^-norm for the solution ■0" as 



liriip 



h — a d — c 



M-l N-1 



j=U k=0 



(3.15) 



for d — 2, and for d = 3 



-i— E E E I'Pjki 

^ j=0 k=0 1=0 



M N 



(3.16) 



For the stability of the time-splitting spectral approximations (A. 7) for 2D and 
(B.l) for 3D, we have the following lemma, which shows that the total density is 
conserved when W{x.) = 0, and resp. decreased when W{x) > 0, in the discretized 
level. 

Lemma 3.1 The time- splitting spectral schemes (A. 7) for 2D and (B.l) for 3D 
GPE with an angular momentum rotation term are unconditionally stable. In fact, 
for any mesh sizes Ax > 0, Ay > and /S.z > 0, and time step size At > 0, 



\\r\\p < wr-'h < w^'Wp = wmp: = 1,2, 

In addition, ifW{'x.) = in (3.1), then we have 



(3.17) 



(3.18) 



Proof: Follows the hne of the analogous results for the hnear and nonhnear Schrodinger 
equations in [6, 9, 10, 7, 12]. 
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4 Numerical results 



In this section, we first test the accuracy of our new numerical method (A. 7) for 
2D and (B.l) for 3D and compare our numerical results with the analytical results 
reviewed in section 2. Then we apply our new numerical method to study vortex 
lattice dynamics in rotating BEC by changing the trapping frequencies and to gen- 
erate a giant vortex by introducing a localized loss term. Our aim is not to find new 
physics phenomena but to demonstrate the efficiency and high resolution of our new 
numerical method 

4.1 Accuracy test 

To test the accuracy of our numerical method in 2D, we take = 100 and Q, — 0.7 
in (1.3). The initial condition in (1.4) is taken as 

7r2 

We take 7^ = 1.0 and = 2.0 in (1.3) and (4.1). Similar example was used in 
[5, 41] for testing numerical accuracy of different numerical methods for rotating 
BEC. The GPE (1.3) is solved on [-8,8] x [-8,8], i.e. we take a = -8, & = 8, 
c = —8 and 0? = 8. Let ip be the exact solution which is obtained numerically by 
using our method with a very fine mesh and small time step, e.g. Aa; = Ay = ^ 
and At = 0.0001, and ipi'^^'^V''^^^ be the numerical solution obtained with the mesh 
size {Ax, Ay) and time step At. 

First we test the spectral accuracy in space by choosing a very small time step 
At — 0.0001, and solving the problem for each fixed (32 with different mesh size 
Aa; = Ay so that the discretization errors in time can be neglected comparing to 
those in space. The errors ||V'(t) — '?/^^^^'^^'^*^(t)||p at t = 0.5 are shown in Table 1 
for different values /?2 and h — Ax — Ay. 





h 


1/2 


1/4 


1/8 




= 20 


2.68E-2 


6.5E-5 


3.41E-10 




= 50 


0.1315 


2.01E-3 


5.91E-8 


P2- 


= 100 


0.4287 


1.94E-2 


8.89E-6 



Table 1: Spatial discretization errors ||'0(t) - '0(^^'^^'^*)(t)||p at t = 0.5 in 2D. 

Next we test the second-order accuracy in time. Table 2 lists the errors at 
t — 0.5 for different values of P2 and time steps At with a fine mesh in space, i.e. 
Ax = Ay = 1/16. 

Similarly, to test the accuracy of our numerical method in 3D, we take (3^ = 100 
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At 


1/40 


1/80 


1/160 


1/320 


1/640 




= 20 


7.86E-4 


1.95E-4 


4.86E-5 


1.21E-5 


3.02E-6 




= 50 


2.95E-3 


7.24E-4 


1.80E-4 


4.49E-5 


1.12E-5 


P2 


= 100 


7.98E-3 


1.98E-3 


4.76E-4 


1.19E-4 


2.96E-5 



Table 2: Temporal discretization errors - V'^^'''^^'^*^^) Ib^ at t = 0.5 in 2D. 

and Q — 0.7 in (1.3). The initial condition in (1.4) is taken as 

774 

We take 7^ = = 7^ = 1.0 in (1.3) and (4.2), and solve the GPE (1.3) in 3D 
on [—8,8] X [—8,8] x [—8,8]. Again let ip be the exact solution which is obtained 
numerically by using our method with a fine mesh and small time step, e.g. Ax = 
Ay = = I and At = 0.0001, and ijji^^AyAzAt) be the numerical solution 
obtained with mesh size {Ax, Ay, Az) and time step At. Table 3 shows the spatial 
discretization errors ||V'(t) - tlj'^^'''^y'^^'^^\t)\\i2 at t = 0.5 with At = 0.0001 for 
different values f]^ and mesh sizes h = Ax = Ay = Az. Table 4 lists the errors at 
t = 0.5 for different values of and time steps At with a fine mesh in space, i.e. 
Ax = Ay = Az = 1/8. 





h 


1 


1/2 


1/4 




= 20 


5.78E-2 


1.27E-3 


2.38E-8 




= 50 


0.1515 


9.50E-3 


2.45E-6 




= 100 


0.3075 


3.88E-2 


7.08E-5 



Table 3: Spatial discretization errors ||V^(t) - i;'''^^'^y'^'^''^\t)\\i2 at t = 0.5 in 3D. 



At 


1/40 


1/80 


1/160 


1/320 


1/640 


/93 = 20 


1.778E-4 


4.435E-5 


1.108E-5 


2.767E-6 


6.897E-7 


As = 50 


6.266E-4 


1.559E-4 


3.892E-5 


9.718E-6 


2.422E-6 


/?3 = 100 


1.63E-3 


1.0379E-1 


1.0069E-1 


2.5111E-5 


6.265E-6 



Table 4: Temporal discretization errors \\ip{t) - '0(^^'^2''^^'^*)(t)||;2 at t = 0.5 in 3D. 

Prom Tables 1-4, we can draw the following conclusions: (i) The method (A. 7) or 
(B.l) is of spectral accuracy in space and second order accuracy in time, (ii) For a 
given fixed mesh size and time step, when f3d is increasing, the errors are increasing 
too. This implies that when the number of atoms in the condensate is increasing. 
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i.e. (5d is increasing, more grid points and smaller time step are needed in practical 
computation in order to achieve a given accuracy. 

Furthermore, Fig. 1 shows time evolution of the normalization N{ip){t), energy 
Ei3,n{ip){t), angular momentum expectation {Lz){t) and condensate widths for the 
above parameters setup in 2D and 3D. 

From Fig. 1, we can see that: (i). The normalization is conserved in both cases 
which confirms (2.1). (ii). The energy is not conserved in the discretized level, but 
the perturbation is very small, e.g. less than 5% (c.f. (a), (c), (e) and (g) of Fig. 
1). (iii). The angular momentum expectation is conserved when = ly = 1 (cf. 
(a) and (e) of Fig. 1) which confirms the analytical result (2.3), and oscillates when 
1 = 7a; 7^ 7y = 2 (c.f. (c) and (g) of Fig. 1). (iv) The condensate widths dx{t), 
Sy{t) and dr{t) are periodic functions when 73. = 7^ = 1 which confirm the analytical 
result (2.6) (c.f. (b) of Fig. 1), and are not periodic functions when 1 — — 2 

(c.f. (d) and (h) of Fig. 1). 

4.2 Dynamics of a vortex lattice in rotating BEC 

In this subsection we numerically study the dynamics of a vortex lattice in rotat- 
ing BEC by changing trap frequencies. This study was motivated by the recent 
experiment [20] in which the frequencies of trapping potential of a stable BEC were 
changed [20]. One of the most striking observation in the experiment is that the 
condensate contains sheet-hke structures rather than individual vortex cores in the 
dynamics by deforming the static trap [20]. By using the hydrodynamic forms of 
the CPE (1.3) in Thomas-Fermi regime, Cozzin et al. [18] tried to study this phe- 
nomena theoretically. But they did not find the sheet-like structures by changing 
the trap frequencies in their theoretical study. Here we study this phenomena by 
directly simulating the GPE (1.3) using our new numerical method. 

We take d = 2, /32 = 100 and Q = 0.99 in (1.3). The initial data ipo{x) in (1.4) is 
chosen as the ground state of (1.3) with d = 2, Q = 0.99, (32 = 100 and 7^ = 7j/ = 1, 
which is computed numerically by the normalized gradient fiow proposed in [13]. In 
the ground state, there are about 61 vortices in the vortex lattice (c.f. Fig. 2). We 
solve the problem on = [-24,24] x [-24,24] with mesh size Ax ^ Ay = 3/16 
and time step At — 0.001. 

First we study free expansion of the quantized vortex lattice. In general, the size 
of a stable vortex lattice in a BEC is too small to visualize it. In experiments, by 
removing the trap, i. e., letting the vortex lattice freely expands, one can obtain 
an enlarged vortex lattice so as to take a photo for it [20]. Of course, they hope 
the vortex structure doesn't change during the free expansion. Thus theoretical 
study of free expansion is very helpful for experiments. We start with the stable 
BEC and remove the trapping at time t — 0, i.e. choosing 7a; — 7?/ — 

in (1.3). 

Similar numerical study was also carried out in [2] by a different numerical method 
with much less number of vortices in the lattice. Fig. 2 shows image plots of the 
density [■0(x, t)|^ at different times for the free expansion of the vortex lattice. From 



10 



14 r 
12 ■ 
10 
8 
6 
4 
2 
O 
2 



(a)-° 





Fig. 1: Time evolution of the normalization N{t) := N{tlj){t), energy Efs^Q{'ilj){t), 
angular momentum expectation {Lz){t) (left column), and condensate widths dxit), 
Sy(t) and 6r(t) (right column). Results in 2D for = 1y = ^ (a&b), and 7^ = 1 
and 7y = 2 (c&d); and results in 3D for 'jx = 1y = Iz = 1 (e&f), and = 1, = 2 
and 7z = 1.5 (g&h). 



the figure, we can see that when the trap is removed at t = 0, the vortex lattice 
will expand with time and the vortex structure as well as the rotational symmetry is 
kept during the expansion. This gives a numerical justification for the free expansion 
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used in BEC experiments. 




Fig. 2: Image plots of the density on [-18,18] x [-18,18] at different 

times t — 0, 0.75, 1.5, 2.0 and 2.75 (from left to right) for the free expansion of a 
quantized vortex lattice. 

Next, we study dynamics of the quantized vortex lattice by changing the trap 
frequencies. We study six different cases: I. 7a; — 1, — 1.5; II. jx = 1)7?/ = 0.75; 
III. 7, = 1.5, 7^ = 1; IV. 7, = 0.75, 7^ = 1; V. 7, = v^, 7^ = V^; VI. 
7a; = v/l-4, -jy = VO.6. Similar numerical study was also carried out in [2] by a 
different method with much less number of vortices in the lattice. Fig. 3 shows 
image plots of the density |t/'(x, t)p at different times for cases 1-11 for changing 
frequencies in y-direction only. Fig. 4 shows similar results for cases III-IV for 
changing frequencies in x-direction only, and Fig. 5 for cases V-VI for changing 
frequencies in both x- and y-directions. 




Fig. 3: Image plots of the density |'(/'(x, t)p on [—18, 18] x [—18, 18] at different times 
t=0, 2.0, 4.0, 6.0, and 8.0 (from left to right) by changing frequency in y-direction 
only from 7y = 1 to: (a) 'jy = 1.5; (b) 'jy = 0.75. 

In Figs. 3-5, initially the condensate is assumed to be in its ground state which 
is a vortex lattice with about 61 vortices. From the numerical results presented 
here, when the trap frequencies are changed at t = 0, we find that: (i) Cases I & II 
correspond to changing trap frequency in y-direction only. The condensate initially 
starts to contract (c.f. (a) of Fig. 3) or expand (c.f. (b) of Fig. 3) in y-direction 
since the trap frequency in y-direction is increasing or decreasing at t = 0. (ii) 
Cases 111 & IV correspond to changing trap frequency in x-direction only. Similar 
results are observed (c.f. (a) and (b) of Fig. 4). (iii) Cases V & VI correspond 
to increasing and decreasing the trapping frequencies in x and y-directions by the 
same value, i.e., e, respectively [18]. The condensate initially starts to contract and 
expand in x and y-directions respectively (c.f. Fig. 5). (iv) We numerically observed 
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Fig. 4: Image plots of the density |'?/'(x, t)p on [—18, 18] x [—18, 18] at different times 
t=0, 2.0, 4.0, 6.0, and 8.0 (from left to right) by changing frequency in x-direction 
only from 7a, = 1 to: (a) 7^ = 1.5; (b) 7^. = 0.75. 




Fig. 5: Image plots of the density t)P on [—18, 18] x [—18, 18] at different times 
by changing trapping frequencies from 7^; = 1 , 7^^ = 1 to 73, = \/l + e, 7^, = y/l.O — e. 
(a) With e = 0.2 for times i = 0, 2, 3, 4, 5, 6, 7 and 8 (from left to right); (b) with 
e = 0.4 for times t = 0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0 and 4.5. 

the remarkable sheet- like vortices in our numerical results (c.f. Figs. 3&5). One 
can compare our numerical results in Figs. 3&5 with the experimental results, e.g. 
Fig. 4 in [20], and find very good qualitative agreement in sheet-like vortex lattice 
formation. Furthermore, we also found that when e is bigger in Cases V&VI, the 
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sheet-like vortices appear earlier. 



4.3 Generation of giant vortex in rotating BEC 

In this subsection we numerically generate a giant vortex in rotating BEC from its 
ground state by introducing a localized loss term [21]. This study was motivated 
by the recent experiment [21] and theoretical study [36] in which the giant vortex 
formation arises as a dynamics effect. We take d — 2, P2 — 100, Vp{x) = and 
Q = 0.99 in (3.1). The initial data V'o(x) in (3.3) is chosen as the ground state 
of (1.3) with d = 2, Q = 0.99, /32 = 100 and 72: = 7y = 1, which is computed 
numerically by the normalized gradient flow proposed in [13]. The localized loss 
term in (3.1) is chosen as a Gaussian function of the form [36] 



W{x, y) = wq exp 



{x - xpf + {y- yof 
rl 



(x,y)eR^ (4.3) 



where Wq, Xq, yo and Tq are constants. We take Wq — 1 and Tq = -\/7/2 in (4.3) 
and solve the problem (3.1)-(3.3) on = [-24,24] x [-24,24] with mesh size 
Ax = Ay = 3/16 and time step At = 0.001. Fig. 6 shows image plots of the density 
|'0(x, t)p at different times for different values of {xo,yo). 

Prom Fig. 6, we can see that the giant vortex lattice is generated due to the 
dynamic effect in a rotating BEC. The center of the giant vortex is the same as 
the center of the localized loss term and the size of the giant vortex depends on 
the values of tq and wq. One can compare our numerical results in Fig. 6 with the 
experimental results, e.g. Fig. 1 in [21], and the theoretical study, e.g. Fig. 1 in 
], and find very good qualitative agreement in giant- vortex formation. 



5 Conclusions 

We have proposed a new time-spitting Fourier pseudospectral method for computing 

dynamics of rotating BEC based on an efficient approximation of CPE with an 
angular momentum rotation term. The new method is explicit, unconditionally 
stable, and of spectral accuracy in space and second order accuracy in time. It is time 
reversible and time transverse invariant in the discretized level, just as the original 
GPE does, and conserves the total density in the discretized level. The efficient and 
accurate numerical method was applied to study dynamics of a quantized vortex 
lattice in rotating BEC. In the future, we plan to extend the idea for constructing 
the new numerical method for one component rotating BEC to multi-component 
rotating BEC and spinor dynamics in a rotational frame, and apply the method to 
study vortex line dynamics in rotating BEC in 3D. 



Appendix A. Discretization in 2D 
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Fig. 6: Image plots of the density on [-12,12] x [-12,12] at different 

times t=0, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75. and 2.0 (from left to right) with wq = 1 
and ro = ^ in (4.3) for generating giant vortices, (a) xq — 0, yo — 0; (b) xq — 1.5, 
yo = 0; (c) xo = 1.5, yo = 1. 

For each fixed y, the operator in the equation (3.11) is in x-direction with con- 
stant coefficients and thus we can discretize it in x-direction by a Fourier pseu- 
dospectral method. Assume 

M/2-1 

'^(x,y,t)^ J2 '$p(y,t) ex.p[iiip(x - a)], (A.l) 

p=-M/2 
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where fip 



2pTT 



b—a 



and ipp{y, t) is the Fourier coefficient for the p-th mode in x-direction. 



Plugging (A.l) into (3.11), noticing the orthogonahty of the Fourier functions, we 
obtain for < p < f - 1 and c < y < d: 



(A.2) 



^jlp + Vtyilpj {t - tn) l/jpiy^tn), tn<t< tn+l- (A.3) 

Similarly, for each fixed x, the operator in the equation (3.12) is in y-direction 
with constant coefficients and thus we can discretize it in y-direction by a Fourier 
pseudospectral method. Assume 



i dtijjpiy, t) = y-Hp + ^yiJ,pj ^p{y, t), tn<t< tn+i- 
The above linear ODE can be integrated in time exactly and we obtain 
V'p(y, ^) = exp ^ '^"^ 



N/2-1 

'^(x,y,t)^ i;q(x,t) exp[iXq{y - c)], 

q=-N/2 



(A.4) 



where = |^ and ■0q(x, t) is the Fourier coefficient for the g-th mode in y-direction. 
Plugging (A.4) into (3.12), noticing the orthogonality of the Fourier functions, we 



obtain for — ^ < q < ^ — 1 and a < x < b: 



i dtiJq{x,t) 



1 



Xl - ilxXq ) lljq{x, t), tn<t < tn+l- 



(A.5) 



Again the above linear ODE can be integrated in time exactly and we obtain 

1 



i>q{x,t) = exp 



^Aj - nxXq ) {t - tn) i>q{x, t„), t^ < t < tn+l- (A.6) 



From time t — tnto t — tn+i, we combine the splitting steps via the standard second 
order Strang splitting [37, 4, 12]: 



Vjk 



M/2-1 

^ ^-^iAt(4+2ny,,l,,)/4 ^^n)^ ^i^xj-a) ^ < J < M, < A; < A^, 

p=-M/2 

N/2-1 

J2 ^-iM{Xl-2nx,X,)/A ^^(l)s) ^iX,{yu-c)^ < k < N, < j < M, 

q=-N/2 " 



^(2) 
jk 



W{xj,yk) = 0, 
, W{xj,yk)>Q, 



N/2-1 

J2 ^-iM{Xl-2Qx,X,)/A (^j3)^ ^iX,{yu-c)^ 0<k<N, < j < M, 

q=-N/2 ^ 

M/2-1 

^ g-iAt(M^+2Qy,^,)/4 ^^(4)^^ ^in,{.x,-a)^ < j < M, < k < N; (A.7) 

p=-M/2 
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where for each fixed k, {'ip^)p {p = —M/2, ■ • ■ , Af/2 — 1) with a an index, the Fourier 
coefficients of the vector = {ipQi^, ipff.^ ' ' "^(m-i)*;)"^' defined as 

- 1 M M 

(V^a = M E V'i^ p = -Y,---,Y-i; (A.8) 

j=0 

similarly, for each fixed j, {'ipf)^ {q = —N/1, • ■ ■ , N/2 — 1), the Fourier coefficients 
of the vector = {ipf^, ip'ji, ■ ■ ■, 'ipf^N-i))'^ ^ defined as 

1 ^-1 M M 

For the algorithm (A. 7) presented in Appendix A, the total memory requirement 
is 0{MN) and the total computational cost per time step is 0{MN\a{MN)). The 
scheme is time reversible when Vr(x) = 0, just as it holds for the GPE (1.3), i.e. the 
scheme is unchanged if we interchange n n + 1 and At <-> —At in (A. 7). Also, 
a main advantage of the numerical method is its time-transverse invariancc when 
W^(x) = 0, just as it holds for the GPE (1.3) itself. If a constant a is added to 
the external potential V , then the discrete wave functions 'ijj'^^'^ obtained from (A. 7) 
get multiplied by the phase factor g-^oC^+i)^*^ which leaves the discrete quadratic 
observable iV'jfc'^^P unchanged. This property does not hold for the finite difference 
scheme [25, 39], the leap-frog spectral method [41] and the efficient discretization 
proposed in [5] for GPE with an angular momentum term. 



Appendix B. Discretization in 3D 

For each fixed y, the operator in the equation (3.13) is in x and z-directions with 
constant coefficients and thus we can discretize it in x and ^-directions by a Fourier 
pseudospectral method. Similarly, for each fixed x, the operator in the equation 
(3.14) is in y and 2;-directions with constant coefficients and thus we can discretize 
it in y and z-directions by a Fourier pseudospectral method. The discretizations of 
(3.13) and (3.14) are similar as those for (3.11) and (3.12) respectively and they are 
omitted here. For simplicity and convenience of the reader, here wc only present 
the algorithm for 3D GPE with an angular momentum rotation term (1.3) with 
< J < M, < A; < A^ and < / < L: 

M/2-1 L/2-1 

^(1) = ^ ^ g-«Ai(2/i2+7j+4f2i/feAip)/8 (^-^ ^iiip{xj-a) ^{■ysizi-e) ^ 

p=-M/2 s=-L/2 

N/2-1 L/2-1 

q=-N/2 s=-L/2 



17 



(3) 
jkl 



-zAtlV{x,,y^,Zi)+Mi^fJl\''] „/,(2) 

^ rjkl 5 

/,(2) 



V', 



^Atwi.^,y,— exp [-i (At\/(,Tj, zi) + 



W{xj,yk,zi) = 0, 
W{xj,yk,zi) > 0, 



iV/2-1 L/2-1 

E E e 

q=-N/2 s=-L/2 

M/2-1 L/2-1 ^-—^ 
p=-M/2 s=-L/2 



where for each fixed (V'^)^^ (-M/2 <p< M/2 - 1, -L/2 < s < L/2 - 1) with 
a an index, the Fourier coefficients of the vector ■0"^^^ (0 < j < M, < Z < L), are 
defined as 



ML ^ 



2 2 ' 2 - 2 

N/1 <q< N/2 - 1, -L/2 <s< L/2-1) with 
a an index, the Fourier coefficients of the vector ■^^^^ {k — 0, ■ ■ ■ , N, I = 0, ■ ■ ■ , L), 

are defined as 



similarly, for each fixed j, {'ipf)^^ 



^ J:Y.^%i e--(--), < g < ^, < . < ^; 



1' NL 



m=0 /=0 



with 7, = 1^ for s = -L/2, L/2-1. 

For the discretization in 3D, the total memory requirement is 0{MNL) and the 
total computational cost per time step is 0{MNL\n{MNL)). Furthermore, the 
discretization is time reversible and time transverse invariant in the discretized level 
when iy(x) = 0. 
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